M)A1  37782 


CR  84.009 


NAVAL  CIVIL  LNCINLLRING  LABORATORY 
Port  iiuencme,  California 


Sponsored  by 

NAVAL  FACILITIES  ENGINEERING  COMMAND 


USER’S  MANUAL  FOR  SAC-3.  A  THREE-DIMENSIONAL  NONLINEAR, 
TIME  DEPENDENT  SOIL  ANALYSIS  CODE  USING  THE  BOUNDING 
SURFACE  PLASTICITY  MODEL 


December  1983 


An  Investigation  Conducted  by 
UNIVERSITY  OF  CALIFORNIA,  DAVIS 


N62583-83-M-T062 


DTIC 


Approved  for  public  release;  distribution  unlimited. 

OTIC  file  copy 


tlCumrv  SIAM*  tic  Avion  OA  TNit  »a0C  Omm  *»*•••*) 


i  AENOAT  DOCUMENTATION  FACE 

READ  INSTRUCTIONS  f 

BEFORE  COMPLETING  FORM  8 

4  TiTcC  feme  SwA4itle> 

User's  Manual  for  SAC-3:  A  Three- 
Dimensional  Nonlinear,  Time  Dependent 
Soil  Analysis  Code  Using  the  Bounding 
Surface  Plasticity  Model 

A  Tv»t  on  M»0«T  A  NENtOO  COVtRCO 

Fipal 

Jan  1983  -  Oct  1983 

A  »C«NO«MING  ONO  NANONT  MUMBCR 

T  >yr.O  «f»i 

Kvran  D.  Mish 

Leonard  R.  Herrmann 

T  cohtk.ct  o*  s*««'  SumRSTm 

N62583-83-M-T062 

University  of  California,  Davis 

IS  MociMliUi.r  msjict  TtH 

YF023.03.01.002 

Naval  Civil  Engineering  Laboratory 

Port  Hueneme,  CA  93043 

*2  «C»ONT  0  A  T  | 

December  1983 

n  NUMACA  ON  NAOIS 

11  MOMI TORINO  a6Cn£v  nan!  A  AOONCIVif  d«f#ar ear  tr#a»  Cmntrmlhmg  Office) 

Naval  Facilities  Engineering  Command 

200  Stovall  Street 

Alexandria,  VA  22332 

If  ICCuNlTV  CV.  All.  ft  *A*s  m»>N) 

Unclassified 

••  OfST*lltUr<ON  STArfMCMT  f«f  **«*  JNp mrtt 


Approved  for  public  release;  distribution  unlimited. 


»T  OtS^AltUTIOM  $▼  ATCmCmT  ft  «*trr«cr  sntv'stf  I*  «f*c*  30.  •!  differs**  trmm  ***•* ) 


»A  »U»*LCMf».TAAv  HQ*€i 


'•  Atv  #0»0S  'CoM'ing*  •*  r»»wn  if  nersaeary  m4  identify  4v  Slock  nwr«i'/ 

Finite  element,  computer  program,  geotechnical  engineering, 
soil  constitutive  law 


20  ABSTRACT  rCMilMt  en  reverse  i(#e  If  necessary  and  identify  Ar  Wsc*  nienAecj 

-The  equations  governing  the  consolidation,  and  the  stress  and 
strains  states  for  soil  structures  are  reviewed  and  their 
historical  development  is  discussed.  Numerical  analysis  con¬ 
cepts  are  used  to  express  these  equations  in  incremental  form. 

A  variational  statement  of  these  incremental  equations  is 
formulated  and  used  in  the  development  of  a  comprehensive 

00. U7J  f O.T.O.  or  .  -ov ••  ’« omol.t.  Unclassified 

SlCuR*Ty  Classification  O*  t«i*  RACt  cNfcen  Date  (mereO) 


finite  element  analysis.  The  concepts  used  in  developing  the 
variational  statement  are  somewhat  different  from  those  used  by 
most  other  investigators  and  appear  to  offer  certain  advantages 
for  inelastic  formulations.  Finally  results  obtained  with  the 
finite  element  analysis  are  compared  to  known  solutions  with 
good  results. 

For  the  convenience  of  the  reader  the  total  report  on  the  ( 
project  is  presented  in  four  parts.  As  noted  above  a  description 
of  the  consolidation  theory  and  certain  theoretical  features  of 
the  finite  element  analysis  are  described  in  the  body  of  the 
main  report  (CR  84.006).  The  second  part  (CR  84.007)  describes 
the  numerical  evaluation  of  the  incremental  form  of  the  bounding 
surface  model.  Finally  "user’s  manuals"  for  the  2-D  and  3-D 
finite  element  programs  are  given  in  two  additional  reports 
(CR  84.008  and  CR  84.009). 


I  JAM  71 


tO'iiom  of  i  NOvi>iSOlfOt.(re 


Unclassified 

l«Cu»'T*  C  <.»>!. »iCAT,o«.  O'  ’“•V  >‘-l  £>.i. 


TABLE  OF  CONTENTS 


L  INTRODUCTION  1 

n.  INPUT  2 

General  Comments  2 

Input  Records  3 

Al.  Title  Record  3 

A2.  Control  Record  3 

A3.  Gravity  Record  4 

A4.  Nonlinear  Analysis  Record  4 

Block  Input  Sections  3 

Bl.  History  Function  Descriptions  5 

B2.  Material  Properties  Array  6 

B3.  Initial  State  Descriptions  8 

B4.  Node  Information  10 

B5.  Surface  Patch  Generation  Information  11 

B6.  Element  Information  12 

B7.  Node  Point  Specifications  13 

B8.  Solution  History  Segment  Data  15 

C.  End  Record  15 

HI.  OUTPUT  16 

IV.  EXPLANATORY  NOTES  REGARDING  THE  INPUT  17 

General  Comments  17 

Section-By-Section  Comments  17 

Al.  Title  Record  18 

A2.  Control  Record  18 

A3.  Gravity  Record  20 

A4.  Nonlinear  Analysis  Record  21 

Block  Input  Sections  22 

Bl.  History  Function  Descriptions  24 

B2.  Material  Properties  25 

B3.  Initial  State  Descriptions  26 

Node,  Patch  and  Element  Data  28 

B4.  Node  Geometry  Information  29 

B5.  Surface  Patch  Generation  Information  30 

B6.  Element  Information  31 

B7.  Node  Point  Specifications  32 

B8.  Solution  History  Segment  Data  35 

C.  End  Record  36 

FIGURES  37 

REFERENCES  45 

EXAMPLE  PROBLEMS 
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INTRODUCTION 


The  finite  element  code  may  be  used  to  analyze  three-dimensional  quasi  - 
static  soil  problems,  including  consolidation  effects.  The  soil  may  be  modeled 
using  either  linear  elasticity  or  the  ’’bounding  surface  plasticity  model  for  cohesive 
soil".  The  program  is  written  in  modular  form  so  that  other  soil  models  can 
be  easily  incorporated.  The  theory  underlying  the  analysis  is  described  fully  in 
the  accompanying  report  [1], 

This  user's  manual  is  divided  into  six  parts:  Introduction,  Input,  Output, 
Explanatory  Comments,  References,  and  Examples.  The  Input  section  gives  in 
outline  form,  the  sequence  of  information  required  to  describe  the  problem  to 
be  analyzed.  Only  the  briefest  notes  of  explanation  are  included  in  this  section. 
Until  an  analyst  becomes  familiar  with  the  program  he  or  she  will  need  to  refer 
to  the  Explanatory  Comments  section  of  the  manual  where  detailed  explanations 
and  examples  are  given.  The  manual  assumes  general  familiarity  with  finite 
element  methods;  novices  in  this  area  are  referred  to  a  standard  text  such  as 
[2]. 


Accession  For  f 
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Special 


Dist 


n.  INPUT 


General  Comments 

Three-dimensional  problems  require  large  amounts  of  input  data,  and  the 
probability  of  errors  due  to  formatted  input  records  and  multi-line  record  lengths 
is  greater  than  for  the  two-dimensional  case.  Therefore,  with  the  exception  of 
the  title  record,  all  input  to  the  three-dimensional  program  is  performed  with 
listed  directed  READ  statements  (i.e.  free-format).  Those  users  who  are  not 
familiar  with  list-directed  input  are  advised  to  study  the  example  input  files 
carfully,  and  to  consult  the  appropriate  FORTRAN-77  documentation  for  the 
computer  being  used.  In  particular,  list -directed  input  records  may  span  several 
lines  (since  input  is  terminated  not  by  line  boundaries,  but  by  exhaustion  of  the 
input  list),  multiple  spaces  between  input  fields  are  ignored,  and  "nuir*  fields  in 
the  input  record  may  leave  the  value  of  the  variable  unchanged  (as  opposed  to 
setting  it  to  zero).  The  safest  way  to  avoid  errors  is  to  include  explicit  values 
for  all  input  variables  that  are  being  read:  this  will  insure  that  the  correct 
problem  is  being  solved. 

Throughout  the  Inppt  section  of  this  manual,  the  following  convention  is 
used:  fields  within  a  record  are  listed  in  order,  and  each  field  is  described  by: 
Type  (I  =  integer,  R  =  real,  L  =  logical),  Name,  and  a  (short)  Description.  Note 
that  (generally)  a  list-directed  read  statement  will  convert  an  integer  value  into 
a  real  value  if  the  corresponding  input  variable  is  real,  but  will  produce  some 
type  of  error  if  logical  data  is  detected  in  a  numeric  (integer  or  real)  field. 


At.  Title  Record  (18A»); 


Any  information  that  is  to  be  printed  as  the  title  of  the  problem. 

A2.  Control  Record: 


Type 

Name 

Description 

L 

IQUIT 

T(rue)  -  stop  analysis  after 

Ffalse)  -  do  not  stop  mesh  Senera,ion 

R 

GRIDV 

= 

grid  generation  parameter  (0.0  -  "isopara¬ 
metric"  grid,  1.0  -  "Laplacian"  grid) 

I 

IFLOW 

=  1 

0  unsaturated  problem  -  no  water  flow 
1  saturated  problem  -  water  flow 

R 

91 

= 

parameter  controlling  numerical  inte¬ 
gration  in  time 

R 

a 

= 

parameter  controlling  "reduced"  integra- 

tion  of  volume  term 

Upper  bounds  on  dimensions  (used  to  establish  dynamic  storage  allocation, 
all  values  except  NCOFMX  >  0): 


NFUNMX 

IFUNSZ 

MATMX 

NCOFMX 

NPTMX 

NELMX 


_>_  no.  of  history  function  specifications, 
Section  B1 

>_  grestest  no.  (M)  of  points  used  to  describe 
a  given  history  function,  Section  B1 

>_  no.  of  materials.  Section  B2 

no.  of  initial  state  descriptions. 
Section  B3 

>_  the  largest  node  number,  Section  B4 


max.  of 


NDSPMX 


number  of  elements 
(Section  B5) 

number  of  patches  used 
to  generate  surfaces  (Section 
B6) 


no.  of  node  point  specifications.  Section 
B7 


V  _VVt*  ' %  ,Vv\  fc‘*  .v  ,*«  ,  *•  ,*»  /  •  ,•  ,  *  -N  .  •  .  •  .  •  . 


Gravity  Record 


Tffie 


Name 


Description 


magnitude  of  acceleration  due  to  gravity 
history  function  associated  with  g 


angles  (in  degrees)  made  by  the  direction 
in  which  gravity  acts  and  the  three 
positive  coordinate  axes 


history  function  associated  with  0  ,  0 

and02g  Xg  yg 


Nonlinear  Analysis  Record: 

(specification  of  desired  iteration  options) 


Tyge 


Name 


NONLIN 


0.0<B<0.5  = 


UMAX 


Description 


[0  linear  problem  )  terminated 

[2}nonhnear  problem  {is  notj  does  not 

occur 

parameter  controlling  Newton-Raphson 
approximation  (0.0  gives  the  tangent 
stiffness  method;  0.5  gives  the  method 
of  successive  approximation) 

maximum  number  of  interations 
permitted  in  any  single  solution 
increment  (default  value  =  5) 


IRPET 


ITFAC 


ERMAX 


0)  ,  ...  (iteration 

(reform  stiffness  matrix  every; 


K-th  iteration* 


T(rue)  variable  acceleration  factor 
applied  to  solution  vector 
F(alse)  components 

places  limits  of  1/FL  >_  (  )  FL  on  the 
acceleration  factor  when  ITFAC  =  T 
(default  value  =  0.3) 

convergence  criterion  for  the  solution 
vector  (default  value  =  0.01) 


*  after  the  2nd 


Bl.  A  record  with  a  1  (integer  one)  in  the  first  field  and  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section),  followed  by  (if  no  history 
function  specifications  are  required,  this  sequence  of  records  is  omitted 
entirely): 

History  Function  Descriptions:  The  following  cards  are  required  for  each 

distinct  function  (History  functions 
numbered  -3,  -2,  -1  and  0  are  explicitly 
defined  in  the  program,  see  Explanatory 
Notes,  and  thus  no  input  is  required): 


1st  Record: 
Type 


Name 

KINDAT 


2nd  record: 
Iffle 


Name 


Description 

must  be  set  to  0  (zero) 

function  number  (>0) 

number  of  points  used  to  define  the 
function 


Description 


function  value 


lR  tm  time  corresponding  to  Fm 

This  record  includes  as  many  fields  as  necessary  to  specify  the  M  pairs 
of  values  (Fm,tm),  m=l-*-M  which  define  the  history  function. 
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B2.  A  record  with  a  2  (integer  two)  in  the  first  field  and  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section),  followed  by: 


r — — — 

distinct 

material: 

1st  Record: 

Type 

Name 

Description 

I 

KINDAT 

must  be  set  to  0  (zero) 

I 

NMAT 

= 

material  number 

I 

ITYP 

= 

( 1  -  isotropic  linear -elastic 
\2  -  anisotropic  linear -elastic 
(  3  -  bounding  surface  plasticity  n<  .1  for 
cohesive  soil 

R 

ps 

soil  density + 

R 

pf 

= 

fluid  density++ 

R 

r 

= 

bulk  modulus  for  fluid  and  soil  particles 
(default  value  =  106) 

effective  soil  permeability  coefficients 


If  the  acceleration  of  gravity  (g  -  Section  A3)  is  taken  as  unity  then  p 
and  p|  are  unit  weights. 

If  it  is  desired  to  use  "excess"  not  total  pore  water  pressure  then  p^  is 
set  equal  to  zero  -  see  explanatory  comments. 
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A  record  with  a  3  (integer  three)  in  the  first  field  and  padded  with 
sufficient  zeros  (see  Explanatory  Comments  section),  followed  by: 


Initial  State  Descriptions:  The  following  information  must  be  supplied 

for  each  non-triviai  initial  state  (this  section 
is  omitted  if  no  information  is  required)* 


Type  Name 


Description 


I  KINDAT  =  must  be  set  to  0  (zero) 


I  IS  NO 


R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 


inital  State  Number 


initial  effective  stress  distribution  for  o 

of  the  form  o  =  a.  +  a~x  +  a,y  +  a, 
x  1  2  r  4 


initial  effective  stress  distribution  for  a 
of  the  form  oy  =  bj  +  b2x  +  b^y  +  b^ 


initial  effective  stress  distribution  for  a 
of  the  form  az  =  Cj  +  c2x  +  c^y  +  c^ 


initial  pore  pressure  distribution  of  the 
form  h  =  dj  +  d2x  +  d^y  +  d^z 


initial  void  ratio  distribution  of  the  form 
e  =  ei  +  e2x  +  e^y  +  e^z 


NN  tsK  NX 


» f  '■  i' ii  *  .■  *  . 


B3  (Continued) 


Type  Name 


R 

R 

R 

R 


Description 


initial  preconsolidation 
distribution  of  the  form  p  =  f 

f3y  ♦  f4z  ° 


The  initial  state  ofo  =  o  =  o  =  h  =  e  =  p  =  0  is  built 
program  as  initial  state  number  0  fzero).  ° 


% 

$ 

g: 

m 
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pressure 
+  f^x  + 


into  the 


/y  iff/. 


B4.  A  record  with  a  4  (integer  four)  in  the  first  field  and  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section),  followed  by: 

Node  Geometry  Information:  As  many  records  as  are  necessary  to  specify 

locations  for  all  nodes  which  are  not  to  be 
generated  with  the  "surface"  and  "interior" 


generation  schemes. 

ll E£ 

Name 

Description 

I 

KINDAT 

= 

must  be  set  to  0  (zero) 

1 

N 

= 

node  point  number 

R 

X 

= 

x  -  coordinate 

R 

Y 

= 

y  -  coordinate 

R 

2 

= 

z  -  coordinate 

I 

R 

INC 

= 

numbering  increment  ^ 

quantities 
l  associated 

D 

xc 1 

spacing  ratio 

with  the 
^  straight  and 
curved  line 

R 

R 

coordinates  of  some 
point  on  the 

generation 
'  options 

YC 

interior  of  the  / 

circular  arc 

R 

ZC 
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A  record  with  a  5  (integer  five)  in  the  first  field  and  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section),  followed  by  (this  section  may 
be  omitted  entirely  if  no  surface  generation  is  needed): 

Surface  Patch  Generation  Information:  As  many  records  as  necessary  to 

specify  the  patches  required  to 
generate  all  surface  nodes:* 


Type 

Name 

Description 

I 

KINDAT 

= 

must  be  set  to  0  (zero) 

I 

N1  i 

the  numbers  of  the  four  nodes  points 

I 

N2  1 

which  describe  the  quadrilateral  or 
triangular**  surface  patch 

I 

N3  ( 

I 

N* 

I 

LIMi 

= 

number  of  additional  patch  layers  in  1- 
direction 

1 

INC1 

= 

numbering  increment  in  1 -direction 

I 

LIM2 

= 

numbers  of  additional  patch  layers  in 
2-direction 

I 

INC2 

= 

numbering  increment  in  2-direction 

this  surface  patch  generation  scheme  is  similar  in  form  to  the  element 
generation  scheme  used  in  the  two-dimensional  program  SAC-2 


for  a  triangular  patch  the  fourth  node  number  is  set  equal  to  the  first 


A  record  with  a  6  (integer  six)  in  the  first  field  and  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section)  followed  by: 

Element  Information:  As  many  records  as  necessary  to  specify  all 

elements  in  the  system* 

Type  Name  Description 

I  KINDAT  =  must  be  set  to  0  (zero) 

I  N1 

I  N2 

I  N3 

1  N4 

I  N5 

I  N6 

I  N7 


N8 

MN 

= 

material  number  (corresponding  to  the 
appropriate  material  description  of 
section  B2) 

ISNO 

S 

initial  state  number  (corresponding  to  the 
appropriate  initial  state  description  of 
section  B3) 

LIM1 

= 

number  of  additional  element  layers  in 
1 -direction 

INC1 

= 

numbering  increment  in  1 -direction 

LIM2 

= 

number  of  additional  element  layers  in 
2 -direction 

INC2 

= 

numbering  increment  in  2-direction 

LIM3 

= 

number  of  additional  element  layers  in 
3-direction 

INC3 

numbering  increment  in  3-direction 

The  numbers  of  the  eight  node  (joints 
which  describe  the  brick  element.  The 
numbering  can  begin  at  any  of  the  nodes 
and  proceed  in  any  direction,  but  must 
occur  in  a  consistent  sequence  (see  figure 

1) .  Various  degenerate  forms  between  a 
brick  and  a  tetrahedron  can  be  obtained 
by  repeating  node  numbers  (see  figure 

2) . 


The  order  of  the  element  records  need  bear  no  relation  to  the  actual 
location  of  the  elements  within  the  body.  The  order  will  determine  the 
assigned  "element  numbers." 


■tfr  . 


a 


I 


Node  Point  Specifications  (continued); 


Iffie 


Name 


Description 


magnitude**  of  the  specified  |d*^“cementl 

for  the  3-coordinate  direction 

history  function  number  for  flow/pressure 

?  [indicates  that  a  known  |water  ! 

P  «pore  water  pressure  > 

is  specified 


magnitude**  of  the  specified  { 


water  flow 
pore  water  pressure 


rotation  angles  (in  degrees)  defining 
the  rotated  coordinate  system  (associated 
with  rotation  option*) 


coefficients  for  distribution  of  applied 
normal  stress***  of  the  form 


an  =  a,  +  a^x  +  a,y  +  a^j 
(associated  with  condensed 
option*) 


pressure 


coefficients  for  distribution  of  water 
source  ***  of  the  form 
q  =  b.  +  b2x  +  b,y  +  b^z 
(associated  with  Condensed  flow 
option*) 


This  option  is  suppressed  by  setting  all  corresponding  values  to  zero. 

In  all  cases  the  actual  value  of  the  prescribed  quantity  is  the  product  of 
the  "magnitude"  and  the  value  of  the  specified  "history  function" 

The  forces/flows  condensed  from  the  pressure/source  distribution  are 
multiplied  by  the  value  of  the  appropriate  history  function  in  fields  IH., 
IH^,  IHj/IH^  (See  Explanatory  comments  section). 
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BS.  A  record  with  an  8  (integer  eight)  in  the  first  field  and  padded  with 
sufficient  zeros  (See  Explanatory  Comments  section)  followed  by: 


Solution  Histc 


it  Data:  One  record  for  each  history  segment 
into  which  the  incremental  analysis  is 
divided:* 


Name 


KINDAT 


NMIS 


TIME 


Description 

must  be  set  to  0  (zero) 

number  of  solution  (time)  increments  into 
which  the  history  segment  is  subdivided 

time  at  the  end  of  the  history  segment 

incrementing  ratio  controlling  the 
timestep  lengths  within  the  history 
segment  (default  value  =  1.0) 


End  Record 

A  record  with  a  9  (integer  nine)  in  the  first  field,  padded  with  sufficient 
zeros  (see  Explanatory  Comments  section). 


The  above  sequence  of  records  Al  ♦  C  are  repeated  for  each  additional  analysis 
in  the  "stack". 


note  that  the  analysis  begins  at  time  tQ  =  0 


in.  OUTPUT 


The  output  from  the  program  consists  of  an  echo  print  of  material 
properties  and  solution  parameters,  the  generated  node  and  element  data, 
messages  for  detected  data  errors,  and  finally  for  each  time  step  the  problem 
solution.  When  data  errors  are  detected,  the  program  aborts  the  job  after  the 
printing  of  the  input  data  and  proceeds  to  the  next  job  in  the  stack  of  data. 

The  printout  of  the  node  point  specifications  includes  any  concentrated 
node  point  forces  (in  x-y-z  coordinates)  resulting  from  specified  surface  pressures 
and  concentrated  flows  resulting  from  specified  surface  sources. 

The  printed  values  of  strains,  stresses,  etc.,  at  any  given  time  step,  are 
the  values  accumulated  to  that  point  in  time  including  initial  values.  The 
stresses  are  effective  stresses  (tension  positive).  The  pore  water  pressure  (units 
of  stress  -  compression  positive)  will  be  either  total  pressure  or  "excess"  pressure 
depending  on  user  preference,  see  Section  B2  in  part  IV. 

The  headings  for  the  solution  output,  are  self  explanatory  with  the  possible 
exception  of  h,  which  denotes  the  pore  water  pressure. 
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IV.  EXPLANATORY  NOTES  REGARDING  THE  INPUT 

General  Comments: 

It  is  the  responsibility  of  the  user  to  maintain  consistent  units.  The  units 
used  to  describe  gravity  (Section  A3),  and  the  material  properties  (Section  B2) 
must  be  consistent  with  those  used  to  describe  the  initial  state  (Section  B3), 
the  geometry  of  the  body  (Section  B4),  and  the  node  point  specifications  (Section 
B6).  The  solution  is  expressed  in  the  same  units  as  the  input. 

Because  the  bandwidth  NBAND  of  the  simultaneous  equations  is  determined 
by  the  numbering  of  the  nodes,  an  optimal  node  numbering  scheme  is  required 
to  minimize  the  computational  cost  of  a  given  finite  element  analysis.  The 
bandwidth  resulting  from  a  given  numbering  scheme  is  computed  in  the  following 
manner: 


i)  Denote  the  span  for  any  two  nodes  of  a  given  element  as  Nj,  where 
N.  is  equal  to  the  absolute  difference  in  the  node  numbers. 

ii)  Denote  the  maximum  value  of  N.  for  a  given  element  j  as  NEj. 

iii)  Considering  all  elements  in  the  system,  denote  the  maximum  value 
of  NEj  as  NEmax. 

iv)  The  bandwidth  is  then  given  by  the  expression  NBAND  =  (3  + 
IFLOW)  *  (NE  +1) 

max 

Since  NE  is  directly  related  to  the  bandwidth  of  the  simultaneous 

ril  aX 

equations,  in  numbering  the  nodes  it  is  this  quantity  that  should  be  minimized. 


Section-by-Section  Comments: 

The  section  numbers  used  below  correspond  to  the  section  numbers  of 
part  0.  INPUT,  thus,  in  order  to  find  information  concerning  the  input  for  B7 
(Node  Point  Specifications)  the  reader  should  refer  to  Section  B7  below.  In 
addition,  within  a  given  section  items  called  out  in  the  input  are  typed  in  bold. 


For  example  input  items  8  and  UMAX  which  are  required  for  input  Section 
A4  are  type  in  bold  where  they  are  discussed  below  in  Section  A4.  The  theory 
underlying  the  analysis  is  only  superficially  treated  here;  for  a  more  complete 
discussion  the  reader  is  referred  to  [1]. 

Al.  Title  Record 

The  title  serves  to  identify  the  particular  problem  under  consideration. 
This  record  must  lie  on  one  72-column  line. 

A2.  Control  Record 

If  a  T(rue)  value  is  specified  for  IQUIT  the  analysis  terminates  after  the 


mesh  has  been  generated  and  printed.  This  option  should  be  used  for  the  first 
run  of  a  large  problem  in  order  to  avoid  wasting  computer  time  analyzing 
incorrect  data.  If  data  for  several  problems  is  contained  in  the  stack,  the 
program  skips  the  time  history  data  for  the  terminated  job  and  proceeds  to  the 
next  problem. 

For  the  precise  meaning  of  the  grid  generation  parameter  GRIDW  the 
reader  is  referred  to  [3]  (GRIDW  =  1.0  -  w,  where  w  is  defined  in  [3]).  In 
general  a  value  of  GRIDW  =  0.0  is  recommended;  for  those  very  rare  cases 
where  this  results  in  a  singular  set  of  equations  for  the  grid  generation  process, 
a  value  of  .05  is  recommended. 

The  code  IFLOW  distinguishes  between  saturated  conditions  where  water 
flow  occurs  (or  a  potential  for  water  flow  exists  -  ideal  undrained  conditions) 
and  unsaturated  conditions. 

When  IFLOW  =  1  the  soil  density  (or  unit  weight  if  the  acceleration  of 
gravity  is  taken  to  be  unity)  p$  specified  in  Section  B2,  refers  only  to  the  soil 
skeleton  (unsaturated  soil).  The  printed  stresses  are  the  "effective  stresses"  and 
must  be  supplemented  by  the  pore  water  pressure  to  obtain  the  total  stresses. 
If  it  is  desired  to  exactly  model  "ideal  undrained  conditions"  (no  movement  of 


water),  the  effective  permeability  of  the  soil  should  be  set  equal  to  zero  (Section 
B2). 

When  IFLOW  =  0  the  soil  density  pg  must  include  the  mass  (or  weight) 
of  any  water  present  in  a  partially  saturated  soil  (the  pore  water  pressure  is 
assumed  to  be  zero  and  water  is  assumed  not  to  flow).  The  printed  stresses 
are  total  stresses. 

Conditions  where  part  of  the  soil  mass  is  unsaturated  and  part  is  saturated 
can  be  modeled  by  specifying  for  the  unsaturated  soil  a  very  small  bulk  modulus 
T  for  the  water  (and  soil  particles  -  Section  B2). 

The  parameter  6j  determines  the  approximation  used  for  the  time 
derivatives  in  the  governing  equations  (see  [1]);  values  between  0.5  (Crank- 
Nicolson)  and  .67  (Galerkin)  are  recommended  [2];  with  the  latter  value  preferred 
when  solution  oscillation  is  a  problem. 

The  parameter  a  determines  the  finite  element  approximation  used  for 
measuring  volume  change  (see  [1]).  When  IFLOW  =  0  a  value  of  0.0  is 
recommended  unless  solution  oscillation  is  a  problem  in  which  case  a  value  of 
.1  may  be  beneficial.  Except  for  nearly  incompressible  linear  elastic  materials, 
when  IFLOW  =  0  a  value  of  1.0  is  usually  preferable. 

All  arrays  in  the  program  whose  dimensions  are  problem  dependent,  are 
dynamically  dimensioned.  The  values  MATMX,  .  .  .  NDSPMX  contain  information 
for  this  purpose.  All  these  quantities,  with  the  exception  of  NCOFMX,  must 
be  greater  than  zero.  The  values  of  MATMX,  etc.  are  upper  bounds  and  thus, 
unless  it  is  desired  to  absolutely  minimize  storage  requirements,  need  not  be 
equal  to  the  actual  number  of  specified  materials,  etc.  The  quantity  NELMX 
must  be  an  upper  bound  both  for  the  number  of  elements  in  the  system  and  the 
number  of  "patches"  used  in  Section  B5  for  surface  generation  purposes;  it  must 
be  an  upper  bound  individually  to  these  quantities  not  their  sum.  When  specifying 
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the  values  of  NELMX  and  NDSPMX  it  must  be  remembered  to  count  those 
elements  (or  patches)  and  node  point  specifications  which  are  included  by  means 
of  the  generation  options. 

In  the  dynamic  dimensioning  of  the  program,  separate  arrays  were  used 
for  integer  and  floating  point  numbers  in  order  to  avoid  difficulties  for  computer 
that  use  different  word  lengths  for  the  two.  The  program  has  been  coded  so 
that  16  bit  integers  may  be  used  if  desired. 

The  dimensions  of  the  program  are  controlled  by  two  quantities  "long" 
and  "longi"  specified  in  "assignment"  statements  at  the  beginning  of  the  program; 
these  quantities  must  satisfy  the  following  inequalities, 
longi  _>  NFUNMX  +  NPTMX  +  9*NELMX  +  (4+IFLOW)*NDSPMX  +  l 
long  _>  31*MATMX  +  2*NFUNMX*IFUNSZ  +  24*(NCOFMX  +  1)  +  [3  +  4*(3  + 
IFLOW)]*NPTMX  +  (6  +  IFLOW)*NDSPMX  +  LONGEQ 
Where  LONGEQ  is  the  space  set  aside  for  solving  the  system  of  equations  by 
means  of  a  block,  constant  bandwidth  equation  solver.  If  the  bandwidth  of  the 
equations  is  denoted  as  NBAND,  then  the  minimum  value  for  LONGEQ  is  NBAND 
♦NBAND  (only  a  single  equation  would  be  contained  in  each  equation  block);  if 
it  is  desired  to  solve  the  equations  entirely  in  core  then  it  must  have  a  value 
_>  NBAND  *  (3  +  IFLOW)  *NPT.  In  general  it  is  recommended  that  LONGEQ 
exceed  the  minimum  by  at  least  30%.  The  calculation  of  the  bandwidth  NBAND 
is  discussed  in  the  general  comments  at  the  beginning  of  this  part  of  the  manual. 

A3.  Gravity  Record 

Gravity  g  can  be  input  either  in  terms  of  the  acceleration  units  appropriate 

2 

to  the  system  of  units  selected  for  the  problem  (32.2  ft/sec  for  English  units) 
or  in  terms  of  multiples  of  the  acceleration  of  gravity  at  sea  level  (i.e.  g  =  1 
for  a  field  structure);  the  corresponding  meanings  of  p$  and  pj  (Section  B3) 
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would  be  mass  densities  in  the  first  case  and  unit  weights  in  the  second.  That 
is,  the  product  pg  must  have  units  of  weight  per  unit  of  volume. 

The  histories  of  the  magnitude  g  and  direction  (0  ,  9  ,  0  )  of  gravity 

x8  yg  zg 

are  specified  by  the  history  function  numbers  (Section  Bl)  IH^  and  IHg.  A 
pre-existing  gravity  loading  of  a  field  deposit  can  be  modeled  by  initializing  the 
stresses  and  pore  water  pressure  (Section  B3)  to  their  proper  values  and  setting 
IH  =  IHq  =  -3.  The  history  of  the  effective  gravity  loading  on  a  centrifuge 
model  during  "spin-up"  can  be  modeled  by  describing  in  Section  Bl  a  history 
function  corresponding  to  the  centrifuge  velocity  history  for  the  test;  in  the 
case  of  a  fixed  bucket  both  g  and  9  would  vary  with  time,  while  for  a  swing-up 

O 

bucket  only  g  would  vary. 

A4.  Nonlinear  Analysis  Record 

For  a  linear  elastic  problem,  NONLIN  and  all  other  input  quantities  for 
the  record  are  set  equal  to  zero.  For  problems  using  the  bounding  surface 
plasticity  model,  NONLIN  is  set  equal  to  1  or  2  depending  on  whether  the 
analysis  should  be  terminated  or  not,  if  convergence  is  not  achieved  in  a  given 
time  step. 

The  factor  8  determines  the  approximation  to  be  used  for  the  Jacobian 
in  the  Newton-Raphson's  iteration  for  the  nonlinear  problem,  for  details  the 
reader  is  referred  to  [4,  5].  It  is  expected  that  a  value  of  0.0  will  in  most 
cases  give  the  best  results. 

The  frequency  of  updating  the  stiffness  matrix  during  the  iteration  process 
is  controlled  by  the  value  of  IRPET  [4];  for  initial  uses  of  the  program  a  value 
of  zero  would  appear  to  be  appropriate.  The  values  of  ITFAC  and  FL  control 
the  use  of  acceleration  factors  applied  to  the  components  of  the  solution  vector 
[4];  for  initial  use  of  the  program  it  is  suggested  that  ITFAC  =  F.  Finally, 


the  default  value  of  .01  for  the  convergence  limit  ERMAX  would  appear  to  be 
adequate  for  most  problems. 

Block  Input  Section 

The  rest  of  the  input  data  for  a  given  analysis  occurs  in  a  block  mode, 
i.e.,  all  the  data  for  history  functions  (for  example)  is  placed  in  one  block,  all 
the  material  properties  in  another,  etc.  Each  block  is  preceeded  by  an  input 
record  with  a  single  integer  in  the  range  1  thru  9  in  the  first  field.  This  integer 
flag  serves  to  instruct  the  program  as  to  the  type  of  data  that  follows  and  is 
characteristic  of  the  finite  element  codes  which  have  been  developed  at  UCD 
and  which  form  the  basis  for  this  three-dimensional  code.  However,  because  a 
list-directed  READ  statement  (used  for  input  to  the  program)  will  ignore  line 
boundaries,  each  record  containing  the  integer  flag  for  a  new  block  of  data 
must  be  padded  with  sufficient  zeros  so  that  the  record  contains  at  least  as 
much  data  as  the  previous  records  (from  the  previous  block).  To  be  specific, 
consider  the  following  case:  The  input  block  for  node  coordinates  is  finished 
and  input  to  the  surface  patch  data  block  is  to  be  read: 

Case  1:  (error) 

0  65  1.75  2.38  5.39  4  1.0  0.0  0.0  0.0  (last  node  coordinate  record) 

5  (integer  flag  for  patch  data) 

0167231245  (1st  record  for  patch  data) 

The  program  would  read  the  5  (indicating  patch  data  follows),  and  in  addition, 
would  continue  until  it  exhausts  the  input  list  (for  node  coordinates),  somewhere 
towards  the  end  of  the  1st  patch  record.  At  this  point,  the  program  and  the 
analyst  would  have  different  opinions  as  to  what  problem  is  being  solved. 

This  difficulty  is  easily  solved  in  one  of  three  ways.  In  the  first  procedure, 
the  input  record  with  the  integer  flag  is  padded  with  sufficient  zero  fields  so 


that  the  input  list  (for  node  coordinates,  in  this  case)  is  exhausted  and  the  next 
read  statement  finds  the  1st  record  for  patch  data: 

Case  2:  (correct  -  method  //l) 

0  65  1.75  2.38  5.39  4  1.0  0.0  0.0  0.0  (node  coordinates) 

500  0  0  000  00  (flag  for  patch) 

016  7  2  312  45  (patch) 

This  method  is  portable,  easy  to  implement,  and  is  used  in  all  the  problems 
in  the  Example  section.  A  second  method  can  be  used  in  environments  where 
a  mechanism  for  terminating  an  input  list  is  available.  For  instance,  a  (/) 
character  will  terminate  the  read  list  on  many  computers: 

Case  3:  (correct  -  method  #2) 

0  65  1.75  2.38  5.39  4  1.0  0.0  0.0  0.0 

5  / 

01  6  7  2  312450 

This  second  method  helps  to  delineate  the  block  structure  of  the  input  ^ile;  a 
check  of  the  appropriate  FORTRAN  -  77  documentation  will  uncover  any  details 
needed  for  a  particular  computer. 

Finally,  the  third  method  is  the  safest  and  may  be  the  easiest  to  use, 
especially  in  environments  where  the  input  data  is  coded  on  cards:  the  flag 
record  consists  of  the  integer  flag  followed  by  28  zero  fields: 

50000000000000000000000000000 
This  method  is  motivated  by  the  fact  that  the  longest  input  record  is 
29  fields  wide  (Node  Point  Specifications).  Therefore,  a  flag  card  with  the  flag 
in  field  one  and  zeros  in  the  next  28  fields  will  always  separate  any  two  data 
blocks.  For  users  accustomed  to  typical  FORTRAN  formatted  input  records, 
this  method  is  recommended. 


Also  note  that  blocks  must  occur  in  the  order  presented  in  the  "Input" 
section  of  this  manual. 

Bl.  History  Function  Descriptions 

The  time  dependence  of  all  input  quantities  (i.e.,  magnitude  and  direction 
of  gravity,  and  specified  node  point  displacements,  flows,  loads  and  pore  water 
pressures)  are  specified  by  means  of  appropriate  "history  functions".  The  program 
has  built-in  four  such  functions  numbered  -3  -►  0,  i.e., 

i)  IH  =  -3  Specifies  a  unit  value  and  a  zero  incremental  value 

for  all  times,  Figure  3a. 

ii)  IH  =  -2  Specifies  a  zero  value  and  a  zero  incremental  value 

for  all  times  -  Figure  3b. 

iii)  IH  =  -1  Specifies  all  incremental  values  equal  to  1.0.  The 

incremental  values  are  taken  to  be  equal  regardless 
of  the  relative  lengths  of  the  time  steps  specified  in 
Section  B7.  The  resulting  history  functions  for  the 
cases  of  equal  and  variable  length  time  steps  are 
illustrated  in  Figure  3c. 

iv)  IH  =  0:  Specifies  a  step-function  at  time  t  =  0;  that  is  a 

quantity  using  this  function  is  applied  entirely  during 
the  first  solution  increment,  Figure  3d. 

In  addition,  the  user  may  describe,  by  means  of  the  input  to  Section  Bl, 
as  many  more  history  functions  (numbered  1-*-)  as  needed;  an  example  of  such 
a  function  is  given  in  Figure  4. 

For  a  particular  history  function,  linear  interpolation  is  used  to  identify 
the  AF  which  corresponds  to  a  given  time  increment  At.  For  any  solution 
times  beyond  the  last  specified  point  t^  the  final  history  segment  is  extended 


indefinitely. 
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When  a  magnitude  V  and  a  history  function  number  IH  are  specified  in 
Section  B7  (or  A3)  for  some  given  external  agent,  then  in  the  solution  interval 
At  an  incremental  value  of  the  quantity  equal  to  V*AF  is  applied,  where  AF 
corresponds  to  history  function  IH. 

B2.  Material  Properties 

The  units  of  the  material  properties  must  be  consistent  with  the  units 
used  to  describe  the  geometry  of  the  body  and  the  magnitudes  of  the  applied 
loads. 

The  material  number  NMAT  serves  as  an  identifier  for  use  in  Section  B6 
to  assign  a  particular  material  description  to  a  group  of  elements.  In  the  current 
version  of  the  program  three  types  of  material  descriptions  are  permitted,  i.e., 
istropic  or  anisotropic,  linear  elastic  and  the  bounding  surface  plasticity  model 
for  cohesive  soils.  Additional  material  models  can  be  easily  added  to  subroutine 
PROPTY  by  extending  the  two  key  "Block  IF"  statements  as  indicated  in  the 
program  by  comment  statements. 

As  noted  previously  (in  Section  A3),  the  units  of  and  p{  must  be 
compatible  with  the  units  seiected  for  gravity  "g". 

Flow  problems  (IFLOW=l)  can  be  expressed  either  in  terms  of  total  or 


excess  pore  water  pressure.  In  the  first  case  p{  must  be  set  equal  to  the  fluid 
density  (or  unit  weight  -  see  previous  paragraph);  in  the  second  case  it  is  set 
equal  to  zero. 

The  quantity  T  can  be  viewed  either  as  the  combined  bulk  modulus  of 
the  soil  particles  and  the  pore  water,  or  as  a  penalty  number  imposing  an 
assumed  incompressibility  condition  for  these  components  [1,  2,  6],  In  the 


absence  of  experimental  evidence,  the  bulk  modulus  for  water  (3.2  x  Hr  psi, 


2.2  x  109  N/m^)  may  be  used  for  f. 
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The  "effective"  permeability  coefficients  k*..  appear  in  Darcy’s  law+  when 
units  of  pressure  (not  head)  are  used  for'  the  pore  water  pressure;  their 
relationships  to  the  permeability  coefficients  commonly  used  by  civil  engineers 
and  those  used  by  physicists  are  discussed  in  [1]. 

For  isotropic,  linear  elasticity  E  and  v  denote  Young's  modulus  and  Poisson's 
ratio  respectively.  The  linear,  anisotropic  elastic  law  is  written  in  the  form: 
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The  meanings  of  the  several  parameters  describing  the  bounding  surface 
model  are  described  in  detail  in  ref  [7-11];  values  for  particular  soils  may  be 
found  in  [8,  11];  a  summary  of  information  is  given  in  Table  1.  In  order  to 
keep  the  input  for  the  model  exactly  as  described  in  [10],  the  parameter  T  is 
retained  even  though  it  is  a  duplication  of  previous  input  (the  second  value  is 
not  used). 

B3.  Initial  State  Descriptions 

The  information  in  this  section  is  used  to  establish  the  initial  state  of 
the  soil.  The  values  of  h  specified  in  this  section  are  used  directly  to  initialize 
the  pore  water  pressure  in  the  elements  and  indirectly  to  initialize  it  for  the 
nodes;  see  Section  B5.  It  is  extremely  important  to  note  that  a x,  a  and  oz 
are  "effective"  stresses  (total  stress  minus  pore  water  pressure).  It  is  assumed 

tin  terms  of  excess  pore  water  pressure  it  has  the  form 

v  -  -(k  *  —  +  k  *  —  +  k  *  — ) ,  etc. 
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that  the  initial  shear  stresses  t  ,  t  ,  and  x  are  zero.  For  linear  elasticity 

^  j  y  a/. 

problems  the  initial  stresses  may  be  taken  to  be  zero  and  then  the  printed 
stresses  are  additions  to  the  initial  state  (i.e.,  superposition  is  valid).  However, 
if  the  bounding  surface  model  is  used,  an  accurate  initiation  of  the  stress  state 
is  extremely  important. 


It  must  be  remembered  that  o 


a  and  a  are  negative  when  compressive, 

y  ^ 


while  h  is  positive  in  compression.  The  units  of  h  are  those  of  stress.  The 


question  of  whether  h  represents  total  or  excess  pore  water  pressure  is  discussed 


in  Section  A3. 


The  initial  states  are  described  by  means  of  simple  linear  equations  in 


the  coordinates  x,  y  and  z  thus  the  coefficients  -»-f^  are  dependent  on  the 
location  selected  for  the  origin  of  the  coordinates. 


The  specifications  of  the  initial  void  ratio  e  and  the  preconsolidation 
pressure  (positive  in  compression)  Pq  are  only  necessary  if  the  bounding  surface 
model  is  used;  they  are  "internal  variables"  for  that  theory  [7]. 

Node,  Patch  and  Element  Data 

This  program  incorporates  several  devices  to  facilitate  mesh  generation. 
When  using  these  procedures  there  is  a  hierarchy  to  the  data  that  must  be 
observed.  First,  the  locations  of  the  nodes  along  all  edges  of  the  body  must 
be  specified  (using  the  straight  or  arc  line  generation  option,  if  desired).  Once 
the  edges  are  specified,  the  surfaces  (whose  boundaries  are  these  edges)  that 
enclose  the  volume  must  be  generated.  These  surfaces  are  collections  of  2-D 
quadrilateral  "patches".  Once  the  surfaces  that  enclose  the  whole  body  are 
generated,  an  extension  of  the  2-D  "interior  node  generation  scheme"  of  Ref.  3 
is  used  to  generate  the  interior  nodes. 

The  analyst  is  urged  to  study  the  example  problems  carefully  to  become 
aquainted  with  the  mesh  generation  schemes.  It  should  be  noted  that  special 


care  must  be  taken  to  define  all_  the  edges  and  ah  the  surfaces  enclosing  the 
body. 

B4.  Node  Geometry  Information 

The  program  incorporates  three  data  generation  routines  to  assist  the  user 
in  defining  the  locations  of  the  system's  node  points:  the  line  or  arc  generation 
option,  the  surface  node  generation  option  and  the  interior  node  generation 
option.  Not  all  numbers  between  1  and  the  maximum  node  number  NPT  need 
correspond  to  actual  nodes  in  the  body.  This  feature  facilitates  the  use  of  the 
available  node  and  element  generation  options.  If  the  location  of  a  node  is 
prescribed  more  than  once  in  the  input  and  the  locations  are  not  in  agreement, 
the  last  description  is  used.  However,  if  in  a  second  or  later  description  the 
node  number  is  entered  as  negative,  then  the  previous  location  is  used.  The 
utility  of  this  option  is  illustrated  later. 

The  straight  line  or  circular  arc  coordinate  generation  option  may  be  used 
whenever  several  sequential  node  points  lie  along  a  straight  line  or  circular  arc 
in  space.  If  such  a  situation  exists,  it  is  necessary  only  to  enter  the  coordinates 
of  the  initial  and  final  points  of  the  sequence  (denoted  by  N'  and  N,  respectively), 
and  the  values  of  INC  and  D.  The  constant  INC  represents  the  difference 
between  any  two  successive  node  numbers  in  the  sequence,  and  D  defines  the 
ratio  of  the  distances  between  any  two  adjacent  pairs  of  points. 

If,  for  a  node  N,  INC  ^  0,  intermediate  node  points  are  generated  along 
a  straight  line  (XC  =  YC  =  ZC  =  0)  or  a  circular  arc  (XC  jt  0  and/or  YC  k  0 
and/or  ZC  £  0)  between  node  N  and  the  point  described  on  the  preceeding  node 
specification  record  N'.  That  is,  the  coordinates  of  the  points  N'  +  INC,  N'  + 
2*INC,  .  .  .,  N  -  INC  are  each  automatically  found.  For  the  case  of  a  circular 
arc  (flagged  by  the  condition  XC  £  0  and/or  YC  ^  0  and/or  ZC  h  0)  it  is  assumed 
to  pass  through  the  end  points  of  the  sequence  N’  and  N,  and  the  additional 
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intermediate  point  with  coordinates  (XC,  YC,  ZC).  This  intermediate  point  need 
not  be  a  node.  The  node  N  for  which  the  specified  non-zero  value  of  INC 
triggers  the  generation  of  the  line  N'  -  N  can  also  serve  as  the  initial  point  of 
a  line  generated  between  it  and  the  point  described  by  the  next  record. 

The  end  points  of  the  sequence  may  be  entered  in  either  order.  For 
example,  the  segments  illustrated  in  Figure  5  could  be  defined  by  specifying  the 
nodes  in  either  the  order  7  +  22  (INC  =  5,  D  =  2.0)  or  the  order  22  +  7  (INC 
=  -5,  D  =  .5).  The  spacing  of  the  intermediate  points  (nodes  12  and  17)  is 
controlled  by  the  spacing  ratio  D.  A  value  of  D  =  1.0  would  result  in  equally 
spaced  nodes. 

B5.  Surface  Patch  Generation  Information 

This  section  is  a  straightforward  generalization  of  a  two-dimensional  mesh 
generation  scheme  [3]  to  surfaces  in  three  dimensions.  If  the  node  point 
coordinate  section  (B4)  completely  specifies  all  the  surface  nodes  (see  the  example 
section,  Terzaghi's  problem),  then  this  block  of  data  can  be  omitted  entirely. 
Under  usual  conditions,  however,  it  will  be  needed. 

Using  this  option,  nodes  on  a  surface  whose  boundary  nodes  have  been 
specified  can  be  automatically  generated  (the  surface  need  not  be  a  plane).  A 
patch  is  a  4-node  surface  element  (which  can  degenerate  into  a  triangle).  A 
surface  is  described  by  an  assembly  of  patches  arranged  in  a  layered  fashion  in 
each  of  the  two  surface  directions,  see  Figure  6. 

The  two  directions  need  not  be  orthogonal,  or  be  aligned  with  the 
coordinate  directions.  The  surface  of  Figure  6  can  be  generated  in  one  of  a 
number  of  ways  (it  is  assumed  that  the  coordinates  of  the  surface  boundary 
nodes  1,  2,  3,  4,  5,  8,  9,  12,  13,  16,  17,  18,  19,  20  have  been  specified): 


Case  1;  (original  patch  defined  by  nodes  1,  2,  6,  5) 
patch  record: 

0  1  2  6  5  2  1  3  4 

Here  (as  shown  in  Figure  6)  Liml  =  2  and  the  1  direction  points  from  node  1 
to  node  4,  Lim2  =  3  and  the  2  direction  points  from  node  1  to  node  17. 
Case  2;  (same  original  patch  as  in  Case  1) 
patch  record: 

0  1  5  6  2  3  4  2  1 


Here,  direction  1  and  direction  2  are  reversed  from  case  1  (and  as  shown  in 
Figure  6).  Other  possibilities  arise  from  choosing  a  different  patch  as  the 
starting  point.  (Note  that  node  numbers  can  be  entered  in  either  a  clockwise 
or  a  counterclockwise  order.) 

For  simple  complete  "box-like”  bodies,  an  entire  face  can  be  generated 
with  just  one  patch  record.  However,  for  more  complicated  bodies,  several 
records  may  be  required  to  specify  all  the  nodes  on  a  given  face.  For  the  face 
pictured  in  Figure  7,  three  records  (corresponding  to  surfaces  A,  B,  and  C)  are 
required. 

The  number  of  patches  generated  by  a  single  input  record  is  (Liml  + 
l)](Lim2  +  1).  The  total  number  of  patches  defining  the  surface  of  the  body 
must  not  be  greater  than  the  dimension  upper  bound  NELMX  (the  same  arrays 
are  used  for  patch  and  element  generation).  Once  all  the  surfaces  enclosing 
the  body  are  specified,  the  interior  node  generation  scheme  automatically  (without 
prompting  by  the  user)  locates  the  coordinates  of  all  the  interior  nodes. 

B6.  Element  Information 

The  material  number  MN  and  the  initial  state  number  ISNO  must  correspond 
to  the  appropriate  descriptions  given  in  Sections  B2  and  B3.  At  the  time  the 
information  for  the  initial  state  ISNO  is  used  to  initalize  h  for  an  element,  it 
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is  also  used  to  initialize  h  for  the  eight  nodes  describing  the  element  in  question. 
If  different  initial  states  are  prescribed  for  two  adjacent  elements  and  they  give 
initial  values  for  h  which  are  not  in  agreement  at  the  common  nodes,  the  values 
obtained  from  the  element  of  higher  number  prevail.  Because  in  practice  h  is 
continuous  such  ambiguous  situations  should  not  often  arise. 

If  the  body  can  be  divided  into  layers  of  elements,  and  if  the  material 
and  the  initial  state  numbers  MN  and  ISNO  are  the  same  for  several  elements 
within  a  layer  and  for  several  layers,  the  node  numbers  of  these  elements  can 
be  simply  established  by  means  of  the  element  data  generation  option.  To 
generate  a  sequence  of  elements  within  a  single  layer,  node  points  are  specified 
for  the  first  element  only,  together  with  appropriate  values  for  LIM1  and  INC1. 
The  generation  of  layers  in  up  to  3  directions  at  a  time  is  a  straightforward 
extension  of  this  concept  (this  is  similar  to  the  2-D  generation  used  in  Section  B5). 
The  total  number  of  elements  generated  by  any  one  record  is  (LIM1  +  1)*(LIM2 
+  1)*(LIM3  +  1). 

Hence,  under  "ideal"  conditions,  the  element  array  for  an  entire  body  can 
be  defined  with  only  a  single  record  in  Section  B6.  Several  of  the  example 
problems  demonstrate  use  of  the  element  generation  option  under  different 
circumstances. 

B7.  Node  Point  Specifications 

Node  specifications  for  displacements  and  loads  may  be  given  in  terms 
of  global  (x-y-z)  components  or  rotated  (xi“x2"x3^  components.  If  0j  =02=0  ^=0.0, 
the  global  system  is  used.  (The  rotation  convention  is  described  later  in  this 
section.) 

For  each  of  the  three  coordinate  directions,  one  may  specify  the  history 
of  either  a  displacement  (IF  =  1)  or  a  load  (IF  =  0)  by  setting  V  equal  to  the 
magnitude  of  the  applied  quantity  and  IH  equal  to  the  appropriate  history  function 


number  of  Section  Bl.  Specified  displacements  and  loads  are  considered  to  be 
positive  when  they  have  the  same  sense  as  the  positive  coordinate  directions. 
In  addition,  either  the  history  of  the  water  flow  Q  or  the  pore  water  pressure 
h  may  be  specified  by  giving  appropriate  values  for  IH^,  IF^  and  V^. 

In  this  section  four  options  are  available  to  aid  the  analyst  in  data 
preparation: 

A.  Generation  of  a  sequence  of  nodes  with  the  same  specifications 

B.  Specification  of  a  rotated  coordinate  system 

C.  Condensation  of  an  applied  normal  stress  to  node  forces 

D.  Condensation  of  a  distributed  water  source  to  node  flows 

Whenever  a  line  or  surface  exists  where  the  node  specifications  are  the 

same  (and  the  node  numbering  is  regular),  option  A  can  be  used  to  generate  all 
the  corresponding  node  specifications  with  one  record.  This  option  can  be  used 
to  generate  such  node  specification  sequences  by  specifying  the  beginning  node 
number  (KK),  the  final  node  numbers  (KKl,  KK2),  and  the  numbering  increments 
in  the  corresponding  directions  (INCI,  INC2).  To  generate  specifications  along 
a  line,  either  KKl  or  KK2  must  be  nonzero  (and  distinct  from  KK).  To  generate 
specifications  over  a  surface,  both  KKl  and  KK2  must  be  nonzero  (and  distinct 
from  KK).  In  either  case  the  appropriate  numbering  increment(s)  INCI  and/or 
INC2  must  be  nonzero. 

Under  normal  circumstances,  the  primary  unknowns  are  the  displacement 
components  in  the  global  (x-y-z)  coordinate  directions.  Forces  and/or 
displacements  can  be  specified  in  a  rotated  coordinate  system  by  specifying  the 
three  rotation  angles  (in  degrees)  necessary  to  transform  the  global  system  to 
the  rotated  one.  A  positive  rotation  is  assumed  to  act  in  a  right-handed  sense 
with  respect  to  the  axis  of  rotation.  The  cumulative  effect  of  these  rotations 
is  the  transformation  from  global  coordinates  (x,y,z)  to  doubly-primed  local 


coordinates  (x",y",z")  according  to  the  covention  shown  in  Figure  9.  Note  that 
these  rotations  are  not  commutative,  i.e.,  the  rotations  must  be  specified  in  the 
order  given. 

An  applied  pressure  distribution  of  the  form  on=aj  +a2x+a3y+a^z  can  be 
automatically  condensed  into  specified  node  forces  by  inputting  the  coefficients 
for  on  such  that  |ajj  +  |a2l  +  ja3 1  +  ja^ j  b  0.  For  such  a  specification  to  be 
valid,  KK,  KK1  and  KK2  must  be  distinct  (otherwise  there  is  no  area  on  which 
the  pressure  acts).  Similarly,  a  distributed  fluid  source  q  =  bj+b2x+b3y+b^z  can 
be  condensed  into  node  flows  by  specifying  coefficients  b^  such  that  |  b ^ |  + 1 b^  | 

+N+KI  *  °* 

The  sign  convention  for  this  condensation  is  that  the  distribution  is  positive 
(i.e.  pressure  acts  inward,  flows  enters  the  body)  if  an  (or  q)  has  a  positive  sign 
and  the  3  vectors  v^,  y2  and  n  form  a  right-handed  system,  where  v.  is  the 
vector  from  node  KK  to  node  KK1,  y2  is  the  vector  from  node  KK  to  KK2, 
and  n  is  an  outward  pointing  normal  from  the  body  (see  Figure  9),  Vj,  y2  and 
n  will  form  a  right-handed  system  if  (Vj  x  y2)  •  n  >  0.  The  condensation 
routine  works  for  flat  or  curved  surfaces,  but  care  should  be  taken  when  using 
it  on  surfaces  with  unusual  curvature  or  with  degenerate  or  non-convex  shapes. 
The  history  functions  for  the  x,y,  and  z  components  of  the  resulting  node  forces 
are  IH^,  IH2  and  IH^,  respectively:  for  condensed  flows,  history  function  IH^  is 
used.  When  the  pressure  option  is  used,  Vj,  V2,  V3,  IFj,  IF2,  1F3,  and  0,,  02 
and  03  are  ignored.  Similarly,  if  the  flow  condensation  option  is  used,  and 
IF^  are  ignored. 

Any  node  may  have  more  than  one  "Node  Point  Specification"  as  long  as 
the  cumulative  effect  of  the  specifications  is  correct.  For  a  given  node,  if  one 
specification  is  a  force  and  the  other  a  displacement  (in  the  same  direction), 
the  displacement  specification  will  prevail.  If  both  specifications  were 


displacements,  the  second  would  prevail.  In  general,  multiple  specifications  at 
corners  are  common.  If  at  a  given  node  it  is  desired  to  use  two  specifications 
where  one  involves  an  applied  pressure  and  the  other  requires  the  use  of  a  local 
coordinate  system,  then  the  pressure  specification  must  occur  first  in  the  data 
block. 

If  IF.  =  NT  =  0  (i  =  1,2, 3, 4)  at  a  node  and  there  is  no  pressure  or  flow 
prescribed,  for  economy,  no  specification  is  needed  and  one  should  not  be  used. 
B8.  Solution  History  Segment  Information 

The  analysis  is  in  general,  time  dependent  due  to  the  consolidation  process 
and  the  history  dependence  of  the  bounding  surface  plasticity  model.  For  a 
non-flow  problem  (IFLOW  =  0)  the  actual  rate  at  which  time  passes  is  not 
important  (because  the  bounding  surface  model  is  rate  independent),  however, 
for  the  purpose  of  modeling  the  history  effects  it  is  still  convenient  to  think 
in  terms  of  time.  For  a  linear  elastic,  non-flow  problem  the  only  role  of  time 

is  to  represent  the  loading  history;  if  only  final  results  are  desired  then  only 

one  time  step  is  required. 

For  convenience,  the  solution  history  is  broken  into  one  or  more  history 
segments.  One  record  is  required  in  Section  B8  for  each  segment.  The  t».ne 
at  the  end  of  a  given  segment  is  denoted  as  TIME;  it  is  assumed  that  the  first 
segment  begins  at  t  =  0.  The  number  of  time  steps  into  which  a  given  segment 
is  to  be  divided  is  prescribed  as  NMIS.  Within  a  given  "history  segment"  the 
ratio  of  two  successive  time  steps  is  equal  to  the  prescribed  spacing  ratio  D; 
a  value  of  1.0  gives  equal  time  steps.  The  role  of  the  time  step  spacing  ratio 

is  analogous  to  the  length  spacing  ratio  used  in  Section  B4  and  illustrated  in 

Figure  5 

The  selection  of  appropriate  time  step  lengths  is  complicated  by  the  fact 
that  two  distinct  processes  are  involved,  i.e.  water  flow  and  soil  plasticity. 


Thus  a  certain  amount  of  experimentation  with  successively  smaller  time  steps 
will  often  be  required.  In  this  process  several  factors  should  be  considered. 
Abrupt  changes  in  step  size  should  be  avoided  (judicious  use  of  the  spacing  ratio 
D  can  facilitate  smooth  transitions  from  small  to  large  time  steps,  etc).  Abrupt 
changes  in  applied  loads  or  displacements  will  cause  large  flow  gradients  and 
require  small  time  steps.  Further  information  concerning  step  size  for  flow 
problems  is  to  be  found  in  [1,2].  A  certain  amount  of  oscillation  of  the  solution 
is  to  be  expected  and  usually  can  be  tolerated.  A  minimum  of  10-20  steps  are 
usually  required  by  the  bounding  surface  plasticity  model  in  proceeding  from  a 
nearly  hydrostatic  stress  state  to  failure  conditions. 

C.  End  Record 

The  function  of  this  card  is  to  signal  the  end  of  the  problem;  the  program 
then  proceeds  to  the  next  stacked  job  (if  any). 
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EXAMPLE  PROBLEMS 


Example  1:  Initial  State  Specification 

The  following  examples  are  intended  to  illustrate  the  use  of  certain  of 
the  input  features  of  the  program.  The  reader  is  referred  to  ref.  [1]  for  a 
comparison  of  results  to  known  solutions. 

This  example  models  an  8-element  cube  with  a  (generic)  initial  state  given 
by  o  =  o  =  o  h  =  e  =  p  =  1.0  +  2.0x  +  3.0y  +  4.0z.  In  figure  10,  the  mesh 

a.  y  o 

is  viewed  from  the  first  octant  (x,y,z  >  0).  The  bottom  surface  of  the  cube 
is  fixed,  and  all  other  nodes  are  free.  Displacements  (output)  are  not  shown, 
since  they  are  all  zero. 

Each  field  where  a  zero  is  intended  as  input  contains  a  zero,  not  a  blank: 
a  blank  (or  several  blanks)  is  a  field  delimiter  in  list-directed  FORTRAN  input. 

A  static  (trivial  in  this  case)  analysis  is  performed  with  one  time  step, 
which  is  an  arbitrary  positive  number  (1.0  in  this  case).  Since  all  the  boundary 
conditions  are  homogenous,  the  history  function  IH  =  -2  is  used  for  the  specified 
displacements:  no  space  is  wasted  by  specifying  zero  forces  on  the  rest  of  the 
boundary. 

Nodes  on  edges  are  specified  first,  followed  by  six  patch  records  to  define 
the  six  faces  of  the  body.  Elements  are  generated  with  one  element  record. 
The  maximum  number  of  "elements"  (bricks  or  patches)  is  twenty-four  (surface 
patches),  not  eight  (brick  elements),  hence  NELMX  =  24.  This  choice  for  NELMX 
is  somewhat  atypical,  due  to  the  coarseness  of  the  mesh:  for  a  more  refined 
mesh,  the  number  of  brick  elements  would  control  the  size  of  NELMX. 

Finally,  note  the  zeros  on  the  record  that  separates  the  material  properties 
block  from  the  initial  state  specification  block.  The  flag  (three)  is  followed  by 
enough  zeros  to  exhaust  the  input  list  for  the  material  properties  read  statement. 
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This  example  models  a  section  of  a  solid  with  a  hole  (tunnel)  in  it.  The 
numbering  of  one  layer  of  elements  is  shown  in  figure  11,  and  the  whole  mesh 
is  viewed  in  figure  12. 

The  surfaces  that  enclose  the  mesh  include  the  hole,  so  the  last  four 
patch  records  define  the  hole's  surface. 

IQUIT  =  TRUE  on  the  second  record,  so  only  the  geometric  (as  opposed 
to  structural)  analysis  is  performed.  This  is  always  a  good  idea  for  the  first 
run  of  a  problem,  especially  when  (as  in  this  case)  the  mesh  is  complicated. 

Note  the  error  message  from  the  program  (excessive  bandwidth)  on  the 
last  output  page.  The  dimensions  of  the  program  as  distributed  are  adequate 
for  testing  purposes,  but  should  be  increased  when  the  user  has  gained  enough 
experience  to  model  large  problems. 
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Figure  12.  Mesh  for  Example  2 
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Example  3:  Terzaghi's  Soil  Column 

This  example  models  a  ten-element  soil  column  with  free  drainage  and 
applied  load  at  the  top  (nodes  41-44)  and  an  impervious  fixed  boundary  at  the 
bottom. 

IFLOW  =  1  in  the  second  record,  since  flows  are  included  in  this  problem. 
Relatively  few  time  steps  are  modelled  (five),  so  the  calculated  dissipation  of 
pore  pressure  is  only  an  approximation  to  Terzaghi's  exact  solution.  Since  excess 
pressure  is  modelled,  the  initial  state  specifications  and  the  densities  are  all 
zero. 

The  applied  load  is  obtained  through  a  pressure  loading  of  the  form 
on  =  10.0.  Note  that  acts  inward,  in  agreement  with  the  sign  convention 
for  the  nodal  force  condensation  option.  (See  Boundary  Specification  section.) 
The  load  is  applied  only  in  the  first  increment,  so  the  corresponding  history 
function  is  IH  =  0.  Note  that  the  pressure  is  applied  before  the  displacement 
specification  for  nodes  (41-44),  as  mentioned  in  the  Boundary  Specification 
section. 

No  patch  data  is  required,  since  the  mesh  is  simple  enough  so  that  all 
exterior  nodes  can  be  defined  without  any  surface  generation. 

Note  that  there  are  44  nodes,  but  that  NDSPMX  =  48  (in  the  second 


record).  This  is  because  there  are  more  nodal  specifications  than  nodes,  since 
nodes  (41-44)  are  specified  twice. 
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Example  4:  Terzaghi's  Problem  Rotated 

This  final  example  demonstrates  the  use  of  the  rotation  convention  for 
boundary  conditions.  The  column  (which  lies  along  the  z-axis  in  example  3)  is 
rotated  into  the  first  quadrant,  where  it  lies  along  the  line  x=y=z.  The  same 
pressure  and  fixed  boundary  condition  specifications  as  those  of  example  3  can 
be  used,  but  the  specifications  for  the  nodes  along  the  sides  of  the  column  must 
be  rotated.  The  rotations  defined  by  0^  =  45.0  degrees,  02  =  35.3  degrees,  and 
0^  arbitrary  (see  figure  8)  are  sufficient  to  rotate  the  column  to  the  desired 
position,  and  are  used  to  specify  the  mixed  force-displacement  condition  for 
nodes  (5-40). 

Although  no  gravity  loadings  were  modelled,  for  illustrative  purposes  the 
angle  in  which  gravity  would  act  was  chosen  to  coincide  with  the  direction  of 


the  column. 
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